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I  3.  ABSTRACT 


Tue  technique  of  second-order  closure  for  turbulence  modeling  is 
used  to  treat  the  growth  and  decay  of  turbulence,  and  the  generation  of 
internal  waves  behind  a  submarine  moving  at  constant  velocity  through 
a  stratified  ocean.  A  computer  mcvicl  based  on  the  solution  of  five 
partial  differential  equatir us  and  approximating  nine  others  is 
generated.  Required  model  coefficients  are  obtained  by  comparison 
with  data  obtained  from  laboratory  flows  and  the  atmospheric  surface 
layer.  The  model  is  verified  by  comparison  with  Naudascher's  wake 
data  in  the  limit  of  v  .nishing  stratif icat ion;  with  Wu's  observations 
in  the  limit  of  starting  with  a  uniformly  mixed  region  in  a  strc  ;gc, 
stratified  fluid;  and  with  Hartman  and  Lewis's  analytical  predictions 
in  the  limit  of  a  small  perturbation  to  the  ambient  linear  density 
gradient  and  vanishing  turbulence.  Model  runs  are  made  demonstrating 
that  the  strength  of  the  internal  waves  generated  by  a  decaying  wake 
increase  with  increasing  initial  potential  energy  (unless  the  initial 
Richardson  number  Ri  «  1  ),  or  increasing  initial  kinetic  energy 
(unless  Ri  »  1  ) ,  or  increasing  ambient  stratification  gradient. 
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SUMMARY 


A  computational  model  has  been  developed  for  t..e  turcui* 

.  submarine  moving  at  a  constant  velocity  through  a  scr tlfi 


'can.  Details  of  the  wake  growth,  collapse  arid  generation  of 


-  rnai  gravity  waves  were  examined  by  an  application  of  the 


:  c*. n  i -order  closure  approach  to  turbulent  flow  developed  at 


> 


...A.p.  over  the  past  few  years.  The  empirical  model  coeffici 
aired  for  basic  Incompressible  flow  were  determined  oy  exumi 
distributions  of  turbulence  in  free  Jets  and  wall  smear  lay 


additional  coefficients  involving  the  influence  of  density 


'atification  were  found  by  comparison  with  observations 
.aspheric  shear  layer. 


The  model  predictions  agree  with  laboratory  observations  in 
.v  limit  of  vanishing  stratification  ana  the  opposite  limit  c’ 
_ ror.giy  stratified  fluid  with  vanishing  turbulence.  Also  in  1 
irr.it  of  a  small  perturbation  to  the  ambient  density  stratlfic: 
::u  vanishing  turbulence,  the  model  results  agree  wi‘h  analytic 
•••edict Ions  of  the  internal  gravity  waves. 

The  influence  of  some  of  the  major  parameters  has  been  de nc 
.  running  a  few  typical  cases.  Increasing  the  Initial  potent l 


rergy  in  the  wake  increases  the  strength  of  collapse  uni  ass 


redact  of  the  Brunt -Valsala  frequency  of  the  fluid  tl.v  s  ’ 
uiracteristic  dimension  of  the  w.-ke  is  much  less  than  tee  cr¬ 
itic  turbulent  velocity  (i.e.,  R1  U,  in  wrier,  case  • :  ®rr 


ffect .  Increasing  tr  Initial  kinetic  energy  incrcast 


collapse  unless  Ri  »  1  ,  in  which  ruse  there  is  no  c:  feci 

c 


creasing  the  Bru.it -Valsala  frequency  of  the-  fluid  •  lw  ^  ?.  ippe 
advance  the  time  and  increase  the  strength  of  collapse. 


■  . 
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model  constants 

drag  coefficient  of  generating  body 

diameter  of  generating  body 

wake  Froude  number  =  r^/U 

gravitational  acceleration 

Brunt -vaisala  frequency  =  [-(g/p)  Sp  /dzj^ 

pressure 

square  root  of  twice  the  turbulent  kinetic  energy 
wake  radius 
initial  wake  radius 

Reynolds  number  =  Ur./v 

1  *• 

Richardson  number  of  turbulence  =  r?N2/q2 

i  m 

model  constant 

velocity  departure  in  free  stream  direction 
normalized  by  U 

Cartesian  velocity  components 

free  stream  uniform  velocity 

horizontal  velocity  normalized  by  U 

model  constant 

vertical  velocity  normalized  by  U 

Cartesian  coordinates 

coordinate  in  free  stream  direction 

coordinate  in  horizontal,  vertical  direction 
normalized  by  initial  wake  radius 

normalized  plotting  coordinate  =  cy,  cz 
(c  =  1.45  in  all  figures  except  Fig.  3  and  6 
where  c  =  2.75) 

iii 
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ir 
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coefl  .cient  of  laminar  diffusion  for  density 
perturbation 

microscale  of  turbulent  dissipation 
mac rosea le  of  turbulent  model 
kinematic  viscosity 
perturbation  pressure  =  p  v  j  gpQdz 
density 

normalized  perturbation  density  =  (p  -  po)/(r1dpQ/dz ) 

ambient  fluid  density 

stream  function  =  /  vdz  =  -  /  wdy 


superscripts 

— —  denotes  time  average 

'  denotes  fluctuation  about  the  mean  value* 


subscript 

m  denotes  maximum  value 
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1 .  INTRODUCTION 

A  submarine  moving  at  constant  speed  through  a  stratified  oceer 
leaves  behind  a  turbulent  wake  containing  potential  and  kinetic 
energy.  The  build-up  of  potential  energy  is  caused  by  the  co-mixing 
of  heavier  density  fluid  above  the  lighter  density  fluid.  The 
kinetic  energy  is  composed  principally  of  the  turbulent  contribution. 
In  the  initial  stages  of  the  wake, the  turbulent  kinetic  energy 
dominates  the  potential  energy  and  the  wake  spreads.  This  spreading 
in  turn  ’ncreases  the  potential  energy,  until  at  some  point  gravita¬ 
tional  forces  are  able  to  overcome  the  turbulence  and  force  a  col laps 
of  the  vertical  spread  of  the  wake. 

The  wake  collapse  has  been  studied  both  experimentally  (Refs,  i-- 
and  theoretically  (Refs.  4-10)  to  gain  a  qualitative  understanding  of 
the  phenomenon.  The  theoretical  studies  have  dealt  primarily  with 
idealized,  inviscid  collapse  of  an  initially  mixed  region,  the  excep¬ 
tion  being  Ko's  analysis  (Ref.  9)  which  attempts  to  include  full 
coupling  between  the  dynamics  of  the  turbulence  and  the  collapse, 
albeit  in  a  simplified,  integral  approach. 

The  present  analysis  uses  the  techniques  of  second-order  clos  .re 
for  turbulence  modeling  (Ref.  11)  to  treat  the  dynamic  coupling 
between  the  turbulence  and  the  stratification  in  some  detail.  We 
are  concerned  here  with  the  growth  and  decay  of  the  turbulence  and 
the  generation  of  internal  waves,  but  not  with  the  px*opagation  of 
these  waves  to  distances  far  outside  the  turbulent  wake. 

The  turbulence  model  is  discussed  in  Section  2.  The  computation.* 
scheme  for  dividing  the  wake  into  regions  governed  by  different  physi 
cal  phenomena  is  discussed  in  Section  3-  Numerical  results  are 
presented  in  Section  4. 
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2.  THE  TURBULENCE  MODEL 


The  time-averaged  equations  of  motion  for  a  Boussinesq  fluid 
whose  motion  consists  of  a  mean  and  a  fluctuating  part  may  be 
written  as 


Whether  the  density  variation  is  caused  by  a  temperature  variation 
or  a  variation  in  compost tlon>  when  the  ambient  density  gradient 
is  constant  the  diffusion  equation  for  the  perturbation  density  P ~  Pc 
may  be  written  as 


D(P-Po)  3  /  dP'fo\  3“jP7 

de - -  3srr 


UJ 
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(3) 


Starting  from  the  equations  of  motion  for  the  fluctuating  compon¬ 
ents  of  velocity  and  density, it  is  possible  to  derive  exact  equa¬ 
tions  for  the  Reynolds  stress  correlation  uTuT  and  the  correla- 

t  J 

tions  involving  the  p'  density  fluctuation.  These  may  be 
written  as 


Du-u. 
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Second-order  turbulence  modeling  Involves  deriving  (or  some¬ 
times  assuming)  relationships  between  the  third-order  correlations 
appearing  in  Eqs.  (4)  -  (6)  and  the  derivatives  of  lower- order 
correlations.  We  choose  to  model  the  third-order  correlations  as 
(Ref.  11) 


Dissipation  terms: 
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Pressure  strain  terms: 
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Diffusion  terms 
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Reasonable  values  for  the  coefficients  a,  b  ,  and  v  '  were 

c 

determined  in  Refs,  11  and  12  as  2.5,  0.125  and  0.3,  respectively, 
by  examining  the  turbulent  distributions  in  free  jets  and  wall  shear 
layers.  The  coefficients  A  and  s  which  influence  the  effect  of 
stratification  were  assigned  the  values  of  0.75  and  1.8, respectively, 
in  Ref.  13  to  obtain  good  agreement  with  the  turbulence  distributions 
observed  in  the  atmospheric  surface  layer.  With  this  choice  of  mode] 
coefficients,  model  predictions  for  the  mean  velocity  gradient,  mean 
temperature  gradient,  Richardson  number,  rms  vertical  velocity 
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J 
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fluctuations,  rms  temperature  fluctuations  and  horizontal  heat  flax 
agreed  favorably  with  surface  layer  experimental  observations  over 
the  complete  range  of  stability  conditions.  The  values  cited  above 


for  the  coefficient  , 
present  analysis. 


a 


A  and  s  will  be  used  in  the 


Equations  (1)  through  (6)  represent  a  formidable  set  for  the 


general  three-dimensional  wake.  To  keep  the  problem  manageable 
for  our  relatively  small  computer  (a  Digital  Scientific  Corp.  META-4), 
we  have  simplified  the  equations  in  certain  regions  of  the  downstream 
development  of  the  wake.  Our  approximation  consists  of  neglecting  the 
convective  and  diffusion  terms  in  the  equations  obtained  by  substituting 


- - - - uyi/aiucu  uy  DUUOlUUl'. 

Eqs.  (7)  through  (16)  into  (4)  through  (6)  to  form  a  set  of  algebraic 


relationships  between  the  turbulent  correlations  and  the  mean  flow 


derivatives,  but  retaining  convection  and  diffusion  of  the  turbulence 
by  carrying  the  equation  for  q2  =  uJuT  #  Equations  (4)  through  (6) 


are  replaced  by 


Note  that  the  function  f  has  been  added  to  Eq.  (18)  to  allow  Eq.  (i 

to  be  added  without  making  the  system  overdetermined.  The  term  is 

added  in  such  a  way  that  the  influence  of  the  convective  and  diffusion 

/ 

terms  in  Eq.  (If)  are  distributed  equally  between  the  three  components 
of  q2  .  Also,  the  diffusion  term  in  Eq.  (17)  is  somewhat  simplified 
from  that  which  would  be  obtained  for  u£u£  from  Eq.  (A). 

The  set  of  equations  (1),  (2),  (3),  ( 17 ) »  ( 18 ) ,  (19)  and  (20)  make 
up  what  we  will  call  the  quasi-equilibrium  (QE)  approximation  for 
second-order  closure.  It  is  closely  related  to  what  Mellor  and  Herring 
(Ref .  14)  refer  to  as  mean  turbulent  snergy  closure.  It  reduces  to  the 
same  number  of  differential  equations  but  retains  more  of  the  coupling 
between  the  components  of  the  Reynolds  stress  through  the  algebraic 
relations,  Eqs.  (18)  -  (20). 

Figure  1  gives  a  comparison  of  the  QE  results  with  that  for  the 
full  set  for  the  case  of  an  unstratified,  axlsymmetric  wake.  Results 
for  the  full  set  are  taken  from  Ref.  12  frcm  an  axlsymmetric  calcula¬ 
tion  while  the  QE  results  are  taken  from  a  3-D,  stratified  calculation 

p 

In  the  limit  of  negligible  sti'atif ication.  The  q  distributions  ag~-  •' 
quite  closely  while  the  mean  velocity  departure  from  free  stream  shows 
aLout  a  \0%  difference  when  it  has  been  allowed  to  decay  an  order  of 
magnitude  from  the  initial  conditions.  The  great  simplification  mad' 
possible  by  the  QE  approximation  appears  to  be  ample  justif ication  for 
this  slight  loss  In  accuracy.  Comparisons  of  the  two  methods  have 
also  been  made  in  calculations  of  the  planetary  boundary  layer  which 
show  relatively  little  difference  between  the  two,  but  these  results 
will  not  be  shown  here. 

The  scale  length  A  in  Eqs.  (17)  -  (20)  should  in  fact  be  solved 

by  its  own  differential  equation,  enabling  us  to  determine  the  A 

variation  across  the  wake  as  the  solution  proceeds  downstream.  The 

2 

equation  for  A  would  oe  similar  to  the  equation  for  q  ,  Eq.  (1/'.' 

(A/e  choose  not  to  do  this,  however,  since  the  scale  equation  has  not, 
yet  become  an  integral  pahi  of  the  second-order  turbulence  model  wt 


are  using,  ins  tend,  we  choose  to  use  a  A  geared  to  the  breadth 
the  q  turbulent  wake  region.  For  our  analysis  here  and  in  tnc 
free  jet  examination  (Ref.  12),  we  assume  that  A  is  proportional 
o  the  distance  iron  the  radius  of  the  maximum  q2  (typically 
r  =  0)  out  to  the  radius  where  the  value  of  q2  has  dropped  to  1/U 
its  maximum  value.  Comparison  with  the  axisymmetrie  flows  In  the 
near-wake  region  yield  a  proportionality  factor  of  0.2. 

The  A  that  results  from  this  algebraic  relationship  now 
expands  or  collapses  depending  upon  the  turbulence  within  the  wake. 
Since, when  collapse  occurs,  the  scale  in  the  vertical  z  direction 
can  be  greatly  distorted  from  the  horizontal  y  scale we  have 
chosen  to  compute  separate  scales  along  the  y  and  z  axes.  These 
resulting  A^'s  and  Az's  then  go  directly  into  their  appropriate 
diffusion  terms  in  the  q‘  equation.  The  QE  approximation  and 
ohe  dissipation  term  -Tbq^/A  use  only  one  scale  length.  Since  w< 
anticipate  that  the  collapse  will  yield  an  ellipse-like  wake  region, 
we  have  chosen  the  equation  for  A  as 


o 


Inus,  when  Ay  »  Az  for  complete  collapse,  A  -*•  2AZ  and  enables 
A  to  track  the  details  of  the  vertical  collapse.  Following  Ref. 
we  also  impose  an  upper  limit  on  the  vertical  scale 
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max 


*  2 
Ri  qd 


_  iL 
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where  Ri  is  a  critical 
surface  layer  to  be  Ri" 


Richardson  number  found  from  the  atmospheri 

0,2b. 


COMPUTATIONAL  SCHEME 


'j 


A  typical  turbulent  wake  in  a  stratified  medium  can  be  divided 
into  regions  governed  by  different  physical  phenomena.  At  the  uft 
end  of  the  body,  the  flow  is  dominated  by  the  details  01  the  body 
shape  and  mode  of  propulsion  (this  problem  will  not  be  discussed 
here).  A  short  distance  downstream  of  the  body  in  the  near  wake, 
the  kinetic  energy  of  the  mean  flow  is  converted  into  turbulent 
kinetic  energy.  As  long  as  the  potential  energy  in  the  near  wak 
is  much  smaller  than  the  kinetic  energy,  which  should  be  true  as 
long  as 
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it  should  bt.  possible  to  1  gnort  the  effect  of  stratification  In  ’ 
region.  Thus  if  there  ar<  no  asymmetries  introduced  by  the  bod;,  , 
near  wake  may  be  treated  as  an  axisymmetr  ic.  flow.  Equations  (2) 

4)  were  integrated  for  unstratified,  axisymmetric  flow  in  Ref.  1 
The  model  demonstrated  the  strong  influence  of  net  mean  momentum  o;. 
the  development  of  the  water,  with  the  model  predictions  verified  5 ,, 
comparisons  with  existing  data  for  a  wake  behind  a  self -propel It  1 
and  for  a  wake  with  significant  mean  momentum. 

The  near-wake  legion  caoulu  remain  valid  as  long  as  the  tunuo 
is  strong  enough  to  dominate  the  stratification,  i.e.,  as  lon^  u 


a 


Although  the  stratification  does  not  influence  the  turbulence 
significantly  in  the  near-wake  region,  the  density  distribution  is 
determined  by  the  turbulence.  Thus,  for  stratified  flow  it  is 
necessary  to  carry  the  additional  Eqs.  (3)>  (5)  and  (o).  However, 
as  long  as  the  initial  conditions  contain  only  the  streaming 
velocity,  it  is  possible  to  reduce  £q.  (2)  to  a  scalar  equation 
for  u  .  Assuming  that  diffusion  in  tne  free  stream  direction  is 
negligible  compared  with  that,  normal  to  the  free  stream,  that  the 
pressure  gradient  in  the  free  stream  direction  is  zero,  and  that 
the  mean  flow  departure  from  the  free  stream  is  small,  we  can  redu- 
Eqs,  (2)  and  (3)  to: 


du 

1  o2,. 

du  'w ' 

du ' 

3x 

w v  u 

dz 

“  3T 

dT 

Pr 

dw'p' 

dT7 

Sr = 

W  v  P 

dz 

'  3T 

where 


I f  we  accept  the  quasi-equilibrium  approximation  for  the  turbuienc 
discussed  in  the  last  section,  then  the  set  is  completed  by  taking 
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together  with  the  algebraic  relations  given  in  Eqs  (8)  -  (,20).  in 
these  equat ions , lengths  have  been  normalized  with  respect  to  the 
Initial  wake  ludius  r^  ;  velocities  with  the  free  stream  velocity. 

U  ;  and  tne  perturbation  density  with  -  in  OhQ/dz  .  The  dimension¬ 
less  parameters  are  defined  as  He  =  U  r  /v  ,  a  Reynolds  number; 

2  —  —  ^ 

Fr  =  [  (  -  gn^  SpQ/^z  )/t  in  2  ,  a  Froude  number;  and  Pr  =  k  /v  ,  the 
Prandtl  number.  The  algebraic  relations  between  the  second-order 
correlations  bring  in  no  additional  parameters  so  the  variables  are 
a  function  only  of  Re,  Pr  and  Fr  plus  the  initial  conditions  on  q  , 
u  ,  and  p  .  The  diffusion  tenms  in  the  turbulent  energy  equation 
have  been  permitted  to  be  more  anise  tropic  in  Eq.  (27)  than  the  form 
in  Eq.  (l/).  This  form  permits  a  larger  influence  of  the  stratifica¬ 
tion,  but  we  have  not  been  able  to  firmly  establish  which  form  is 
more  accurate  yet. 

Equations  (2$)  to  (27)  cease  to  be  a  valid  approximation  for 
the  wake  when  the  poten*  Ial  energy  due  to  the  perturbed  density 
Increases  to  the  point  where  it  is  the  same  order  of  magnitude  as 
the  kinetic  energy  of  the  turbulence,  i  .e . ,  when  Ri  >  0.1  .  Then 
the  effect  of  the  induced  vertical  velocity  on  the  convection  terms 
can  no  longer  be  neglected.  To  continue  the  wake  calculations  Into 
the  collapse  region,  it  is  then  necessary  to  include  the  two  transverse 
momentum  equations  from  Eq.  (2). 

In  the  case  of  a  momentumless  wa>ke,  the  mean  velocity  departure 
from  free  stream  decays  much  more  rapidly  then  the  rms  turbulent 
velocity  so  that  with  density’  stiatif ieation  typical  of  the  ocean 
^ Brunt -Vaisala  period  of  the  order  of  one-hall’  hour)  the  mean  velocity 
departure  becomes  negligible  before  Rl  reaches  0.1.  Therefore,  for 
many  cases  of  Interest  we  can  set  u  -  c  and  drop  the  axial  momentum 
equation  In  the  collapse  region.  The  set  of  equations  to  be  solved 
In  this  region  within  the  quasi-equil  ibrlurn  approximation  now  becomes 


OV  dw 


dF  dy  33T[3y  4  ^Fj 


The  set  is  completed  with  the  algebraic  relationships  for  the  correla 
tions,  Eqs.  (18)  to  (kO).  Note  that  rather  than  work  directly  with 


the  continuity  equations,  Eq.  ( 1 ) ,  we  have  taken  the  divergence  of 
the  momentum  equal  !or.  to  obt  ain  a  Poisson  equation  for  the  perturba¬ 


tion  pressure. 


d/. .  The  last  term  In  Eq.  (32)  has  been 


retained  because  previous  researchers  found  1 t  desirable  tc  correct 


for  the  fact  that,  numerically,  Eq.  (i)  cannot  be  satisfied  exactly. 


In  summarizing,  trie  d  l vision  of  the  stratified  wake  development 
into  a  growth  Phase  (!)  gov-  ;  ed  by  Eqs.  ( .-*5)  -  (27)  and  a  collapse 


Phase  (II)  governed  by  Eqs. 


fc  car  s_j  that  (I)  should 


hold  as  long  as  Ri  <  >j.  l.  whl'i-  ;  l  »  sr.o  *ri  no l u  when  |u  l/q  <  0.1 

u;  in  - 


_ 


3-5 


In  the  overlap  region,  which  should  occur  as  long  as  Ri  immediately 
behind  the  body  is  sufficiently  small,  either  set  of  equations  may 
be  useo.  If  there  is  no  overlap,  such  as  may  occur  in  the  wake 
behind  a  low  drag  body  moving  slowly  through  a  steep  stratification, 
then  our  present  computational  scheme  is  not  valid.  Future  modifica¬ 
tions  to  our  computer  model  will  provide  for  the  simultaneous  solution 
of  all  three  momentum  equations  so  that  this  more  general  case  can  be 
analyzed.  Figure  2  gives  a  comparison  between  tne  results  of  a  Phase 
(I)  and  (II)  computation  for  a  case  with  a  large  overlap  region.  No 
significant  differe  ices  are  observed  in  the  two  solutions  until  after 
Ri  =  0.07  at  .04  B.V.  (B:  unt-Vaisa.la  periods)  after  wake  generation. 

Of  course.  Phase  (I)  represents  a  much  more  efficient  computational 
scheme  in  the  overlap  region. 

We  have  all  of  the  information  needed  to  solve  the  equations  of 
interest  for  our  Phase  (I)  or  (lb'  regions.  We  begin  our  calcula¬ 
tions  In  the  near-wake  region  by  assuming  initial  conditions  fcr  the 
turbulence,  the  perturbation  density  and  mean  u  profile.  The 
Phase  (I)  region  is  permitted  to  build  to  the  overlap  region,  at  which 
time  we  introduce  the  1  trmai  velocities  v  and  w  and  the  pressure 
v  to  proceed  through  the  collapse  Phase  (ll).  To  reduce  redundant 
calculations,  we  consider  only  a  quarter-plane  of  the  wake  cross- 
section,  thus  requiring  us  to  impose  appropriate  symmetry  conditions 
along  both  axes. 

Since  our  primary  purpose  is  to  calculate  the  wake  dynamics  rather 
than  the  full  wave  pattern  between  the  wake  and  the  surface,  we  impose 
an  artificial  porous  liner  a  few  radii  from  the  turbulent  wake;  this 
liner  serves  to  absorb  any  waves  emitted  from  the  wake  without  reflect¬ 
ing  them  back  to  the  wake.  In  this  way  we  are  justified  in  requiring 
all  the  variables  -*■  C  as  r  -►  ^  .  This  absorption  of  the  waves  is 
implemented  by  adding  ih-3  damping  terms  -  kw  and  -  kv  to  the  right- 
hand  sides  of  Eqs.  i2B)  and  (  y },  respectively,  and  the  term  -  v  dk/dy 
-  w  dk/dz  to  the  prsssur-  nq.  (32).  The  damping  factor  k  is  zero 
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out  to  a  few  wak  radii  and  then  increases  smoothly  to  absorb  any  waves 
A',  the  wake  dimension  grows,  the  position  of  the  liner  also  moves  out 
at  a  proportional  rate. 

Figure  3  is  a  plot  of  density  contours  for  two  different  positions 
of  the  absorbing  liner.  Although  the  contours  are  quite  different 
outside  the  inner  liner,  there  is  little  difference  in  the  region  of 
interest  around  the  turbulent  portion  of  the  wake  itself.  A  detailed 
comparison  of  these  two  cases  shows  that  the  average  numerical  differ¬ 
ence  at  the  points  inside  the  inner  liner  is  less  than  of  the 
maximum  local  value. 

The  differential  equations  discussed  above  are  solved  numerically 
by  replacing  them  witn  i'orward-time-centered-space  finite  difference 
equations  (Rei  .  1m  .Since  he  equat  ions  involve  two  space  directions 
y  and  z  and  a  t! rue-1  ike  direction  x  ,  we  choose  to  solve  them  by 
the  Alternating  Direction  Implicit  (ADI)  method  (Ref.  18).  The 
implicit  technique  iermHs  s' able  solution  development  with  a 
contrcllec  increase  in  Ax  per  step.  The  mesh  In  the  y  and  z 
directions  is  variable,  w1  Ih  mesh  spacing  Ay  ana  Az  minimizing 
curvature  changes  from  point  to  point. 

When  needed,  the  ioisscn  equation  (Eq.  (32))  for*  the  perturbation 
pressure  t r  is  solved  by  adding  a  time -like  term  -  bv/bt  to  the 
left-hand  side  of  Eq.  (32)  and  iterating  by  the  ADI  method  to  a  steady 
state  solution  (so  that  dv/di  =  0)  for  every  step  in  q^  ,  p  ,  v 
and  w  .  For  the  I  irst,  two  iteration  steps.  At  is  given  a  value 
proportional  to  the  urea  of  tne  wake  (At  =  yjnax  ztnax/^0)  and  allowed 
to  increase  or  decrease  thereafter  at  a  rate  determined  internally  in 
such  a  way  as  to  try  to  maintain  a  change  of  =  10#  with  each  iteration. 

In  practise  our  solution  scheme  takes  a  step  Ax  in  the  main 
variables  (q^  ,  ,  ,  md  u  ,  or  v  and.  w  as  the  case  dictates), 
then  the  pressure  is  iterated  for  these  new  values.  We  find  that 
our  procedure  works  very  veil  for  Phase  (I)  (where  no  pressure 
iteration  is  required),  but  the  Phase  (II)  pressure  solution  has 


t 


generation  for  twc  positions  of  the  porous  liner.  Oi 
sids  olie  indicated  dashed  line  the  liner  resistance 

increases  from  zero  exponentially.  -  contours  for 

the  case  of  the  larger  liner,  (initial  conditions  - 


—  \ - — *  -  - - —  AU  -IWiiU  — 

Fr  =  1.15,  qm=ujjj  =  v  =  w  =  0,  Pm  =  1 . )  fm  =  .28c 
v;ith  contour  for  P/Pm  -  +.5  denoted  by  -1-;  -.5  bj 
+*1  ”3”,  “•!  by  -4-;  +:01  by  -5-;  and  -.01  b 

-0-. 


1  '  f  m  ^  *■  m  _  *  ^ 

denoted  by  -1-;  -.5  b; 


caused  the  bulk,  of  our  numerical  problems.  It  seems  that  our  ADI 
method  is  not  an  extremely  viable  one  when  steady  state  is  sought 
for  t t  at  ever.v  step  Ax  .  Rather,  we  found  that  stability  and 
convergence  were  better  controlled  by  limiting  the  number  of  pressure 
iterations  at  each  step;  a  maximum  of  five  iterations  were  permitted. 
This  serves  to  keep  the  pressure  close  to  an  accurate  value,  and  in 
turn  keep  v  and  w  ,  tnus  p  ,  accurate  also.  We  monitor  our 
accuracy  by  trying  to  keep  the  relative  error 


J  (Aw)2dy  dz 

r  2 

j  F  ay  dz 


i 


2 


<  0.001 


where  Att  is  the  change  in  it  in  the  step  At  ,  and  F  is  the 
value  of  the  right-hand  side  of  Eq.  (32).  F  has  also  been  slightly 
modified  by  realizing  that  the  term  of  most  importance  -  l/'Fr  dp/dz 
can  be  estimated  for  the  next  step  size  Ax  and  the  pressure  estimated 
ahead  accordingly.  Since  the  numerical  technique  is  not  exactly 
divergence  free,  the  extra  term  on  the  right  side  of  Eq.  (32)  in  the 
form  dV-V/dx  is  retained  to  make  V*V  =  0  at  the  future  Ax  main 
step. 

A  typical  run  in  Phase  (I)  takes  about  3  hours  on  the  Digital 
Scientific  METa-4  (12  min.  on  CDC  6600),  while  Phase  (II)  takes 
about  12  hours  on  the  META-4  (30  min. on  CDC  6600)  of  which  about 
50$  of  the  time  is  involved  In  the  pressure  iteration  loop. 
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4.  NUMERICAL  RESULTS 

Considerable  confidence  in  our  turbulence  model  was  gained  by 
the  comparisons  ith  che  unstratified  wakes  in  Ref.  12  and  with 
the  stratified  atmospheric  shear  layer  in  Ref.  13.  To  gain  some 
confidence  in  our  computer  model's  capability  to  predict  the 
generation  of  internal  waves,  we  made  a  check,  with  Hartman  and 
Lewis's  linear,  inviscid  solution  (Ref.  3).  Within  our  self- 
imposed  limit  of  40  x  40  mesh  points,  we  were  unable  to  maintain 
good  agreement  in  the  immediate  neighborhood  of  the  discontinuity 
at  the  edge  of  the  homogeneously  mixed  region.  However,  results 
within  the  center  of  the  mixed  region  agree  with  the  exact  solution 
to  times  greater  than  one  B.V.  Figure  4  is  a  comparison  of  the  model 
results  with  the  linear  analytical  predictions  for  radii  less  than 
80$  of  the  radius  at  which  the  discontinuity  occurs.  The  bars  on 
the  model  results  give  the  r.m.s.  deviation  of  the  value  at 
individual  points  about  their  mean. 

Another  cheek  on  the  model  is  given  in  Fig.  3  by  comparing  with 
the  experimental  results  if  Wu  (Ref.  3).  WU  used  a  mechanical  mixer 
to  produce  a  cylindrical  region  of  homogeneously  mixed  fluid 
surrounded  by  a  strongly  stratified  fluid.  His  observations  of 
the  width  of  the  mixed  region  as  a  function  of  time  after  release 
are  plotted  as  given  in  Ref.  10.  The  model  prediction  for  growth 
of  the  wake  width  for  a  case  with  lai  ge  initial  Ri  and  an  initial 
homogeneously  mixed  wake  is  given  as  the  solid  curve.  The  agreement 
is  similar  to  that  obtained  by  the  numerical  models  of  Refs.  6  and  10. 
The  streamline  patterns  predicted  for  this  two-dimensional  flow  at 
times  of  0.5,  1.0  and  1.5  B.V,  after  release  are  presented  in  Fig.  6. 
At  earlier  times  there  is  only  one  vortex  in  the  quadrant  but  other 
modes  appear  as  time  progresses.  At  t  =  1.0  B.V.  there  are  three 
distinct  vortices,  while  at  1.5  B.V.  two  o'f  these  of  the  same  sign 
have  joined  and  a  fourth  is  forming  along  the  vertical  axis. 


Figure  4.  Perturbation  density 
and  velocity  variations  as  a 
function  of  time  for  an  inviscid, 
small  amplitude  linear  perturbation 
in  density  within  a  cylinder  of 
fluid.  e(=  10“5)  defined  by 
initial  condition  on  p  .  -  ana¬ 

lytic  solution  given  by  Hartman  and 
Lewis  (Ref.  5),  numerical  model 
predictions  with  r.m.s.  scatter  in 
the  inner  Q0%  of  the  radius  of  the 
cylinder  as  indicated. 
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Figure  5. 
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tlons,  V  experimental  observations  of  Wu  (Ref. 
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For  our  case  of  a  momenturnless  wak.e  the  principal  parameters  are 
the  Froude  nurrber  Fr  =.  [  - (  g/ p )  dpQ/dz  ]  2U/D  ,  the  initial  level  of 
turbulence  q  /U  .  the  initial  density  profile  in  the  wake,  and  the 
Reynolds  number  of  the  turbulence  (  q  A/v)^n^t.  For  a  submarine  wake 
the  most  important  single  parameter  is  a  combination  of  the  first  two, 
a  Richardson  number  for  the  stratified  turbulent  wake  (Eq.  (23)).  As 
seen  in  the  previous  section  when  Ri  «  1  right  behind  the  body,  the 
wake  will  grow  for  some  distance  behind  the  body  before  the  wake 
collapses.  Then  the  initial  density  distribution  is  relatively 
unimportant  since  the  turbulent  entrainment  of  fluid  in  the  growth 
region  will  lead  to  a  density  distribution  at  the  point  where  collapse 
begins  which  is  independent  of  the  initial  distribution.  The  Reynolds 
number  is  not  influential  as  long  as  it  is  large  enough  so  that  qA/v 
remains  much  larger  than  one  through  the  initial  collapse.  The  only 
parameter  we  will  consider  in  detail  herein  is  Ri. 

If  Naudascher's  (Ref.  15)  values  for  turbulent  intensities  behind 

a  self-propelled  body  are  used  as  a  starting  point  (although  Gran  (Ref.  lb ) 

indicates  that  his  values  are  about  a  factor  of  two  higher  than  would  be 

obtained  behind  a  streamlined  body  driven  by  an  efficient  propeller),  then 

a  body  with  a  radius  of  five  meters  moving  through  a  stratified  fluid 
2  -5  -2 

with  N  =  10  vsec  at  a  velocity  of  five  meters/sec  would  have 
_2 

.Ri  2  10  at  a  station  20  diameters  downstream  of  the  body.  Figure  7 
is  a  plot  of  some  of  the  wake  variables  as  a  function  of  distance  down¬ 
stream  of  the  body  or, equivalently, time  after  generation  for  a  run  in 
which  we  start  with  a  homogeneously  mixed  density  and  Naudascher's 
values  of  turbulence  at  x/D  -  20  .  Up  to  x/D  s  200  where  Ri  ~  1, 
the  results  for  q  are  the  same  as  obtained  in  Ref.  12  for  an  unstrati- 
fled  wake;  the  density  perturbation  is  being  reduced  as  the  wake  entrains 
more  fluid;  and  v  and  w  are  much  less  than  q  .  At  larger  distances 
Jownstream  there  is  a  significant  difference  in  the  shape  of  the  wake 
although  the  maximum  c  shows  little  difference.  Although  v  and  w 
build  to  maximum  velocities  at  approximately  0.5  B.V.,  they  remain 
smaller  than  q  .  Figures  6,  9  and  10  show  the  contours  of  p  and  q" 


LOG  (MAX  VALUE} 


Figure  7 
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LOG  (x/D -20) 


Maximum  values  of  q  ,  p  ,  v  and  w  as  a  function  of 
time  after  generation.  (Initial  conditions  -  Fr  -  ?28 
=  *00^56,  v  =  w  =  0  pm  =  1  . )  Equivalent  distance 
behind  a  Naudascher  body  Is  also  indicated. 
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Figure  9.  Contours  of  constant  p  at  different 
conditions  as  Fig.  7.  Contours  for 
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and  the  streamline  patterns  of  the  radiating  waves  at  0.0,  0.5>  1.0 
and  1.5  B.V.  periods  after  t lie  wake  generation. 

When  the  stratification  is  increased  the  picture  changes  as  seen 
in  Figs.  11  to  14.  Stratif icaticn  oecomes  important  a  shorter 
distance  behind  the  body,  before  the  wake  has  grown  to  as  large  a 
radius  as  previously,  and  the  strength  of  the  streamline  patterns 
produced  by  the  collapse  is  increased.  As  a  function  of  x/D  ,  q 
and  p  are  the  same  for  both  runs  until  x/D  -  67  (~  0.1  B.V.  on 

A 

Fig.  11)  at  which  point  the  stronger  stratification  forces  p  and 
q  to  be  decreased  at  a  faster  rate  as  v  and  w  increase.  In  the 
region  of  the  principal  collapse  (0.1  to  1  B.V.)  v  and  w  scale 
fairly  closely  with  N  ,  as  seen  in  Figs.  7  and  11,  and  the  stream¬ 
line  patterns  are  quite  similar  as  seen  in  Figs.  10  and  14. 

When  the  initial  Richardson  number  is  decreased  by  increasing 
the  intensity  of  the  turbulence,  the  results  shown  in  Figs.  15  to  1 8 
are  obtained.  The  Influence  of  stratification  is  delayed  until  the 
potential  energy  can  build  to  a  larger  vtlue.  Consequently,  a 
stronger  streamline  pattern  is  generated  by  the  collapse  although 
its  occurrence  is  delayed  further  downstream. 


The  maximum  potential  energy  per  unit  length  of  the  wake  occurs 
when  the  local  Richardson  number  is  of  order  pne.  Figure  19  shows 
how  this  maximum  value  varies  as  a  function  of  the  initial  Richardson 
number  over  the  range  of  initial  Ricnardson  numbers  expected  to  be  of 
interest  for  a  submarine  (IC~H  <  Hi  <  i)  when  the  initial  distribu¬ 
tions  of  density  and  turoulenee  are  held  constant.  Sensitivity  to 
variations  in  trie  diet :i: at ion  of  the  initial  turbulence  have  not 
been  computed  cut  It  I  ex.  cted  tha‘  these  would  have  little  influ¬ 
ence  on  the  aependervj'  anew:.  In  Fig.  19  as  long  as  the  initial 
kinetic  energy  was  fi... -a.  f  t  results  art  sensitive  to  variations 
in  the  Initial  ratio  o:  potential  to  kinetic  energy.  Two  curves 


obtained  by  starting  xltn  u  if’ t.- re¬ 


nt  Initial  density  profiles  are 


Included  on  this  curve  to  demonstrate  in  Is  variation.  As  xong  as 


RiQ  is  small  enougn,  the  initial  density  a  istrj  nut  ion  is  unimportant 
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Maximum  values  of  q  ,  p  ,  v  and  w  as  a  function 
of  time  after  generation.  Initial  conditions  are 
the  same  as  Fig.  7  except  Fr  is  decreased  by  v/Io. 


Figure  12.  Contours  of  constant  turbulent  kinetic  energy  at  different  times  after  the 
generation  for  the  same  conditions  as  Fig.  11.  Contour  notation  per  Fig. 


Figure  13.  Contours  of  constant  p  at  different  times  after  generation  for  the  same 
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Figure  14.  Normalized  streamline  patterns  at  different  times 

after  generation  for  the  same  conditions  as  Fig.  11 
Contour  notation  per  Fig.  3. 
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Figure  1 6.  Contours  of  constant  of  at  different  times  after  generation  for  the 
same  condition  as  Fig.  15.  Contour  notation  per  Fist. 
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since  the  wake  will  continue  to  grow  and  develop  its  appropriate 
density  distribution  prior  to  collapse.  When  RiQ  >  1,  collapse 
is  imminent  and  the  initial  density  distribution  is  a  dominant 
factor.  Since  in  this  case  the  initial  potential  energy  is  the 
maximum  value  of  potential  energy  it  is  independent  of  RiQ  . 

The  limiting  asymptote  for  zero  initial  potential  energy  varies 
approximately  with  the  -0,6  power  of  the  initial  Richardson  number. 

A  summary  curve  showing  the  distortion  of  the  wake  as  a  function 

of  time  after  generation  is  presented  in  Fig,  0  for  several  of  the 

points  on  Fig.  19.  The  experimental  variation  as  reported  by 
Schooley  and  Stewart  (Ref.  1)  is  also  inclidf-d.  A  detailed  match 
of  Schooley  -  Stewart's  run  has  not  been  made  yet  since  their  initial 

Ri  «  0.1  ,  where  it  is  necessary  to  also  know  the  initial  density 

distribution  in  the  wake.  It  appears  that  a  reasonable  match  could 
be  made  since  their  results  fall  between  trie  Ri  =  0.04o5  and  0.465 
model  results. 
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We  also  recommend  th.t  'experimental  laboratory  data  be  chosen 
by  the  contract  monitor  so  that  detailed  comparisons  between  our 
model  predictions  and  observations  can  b<  made  for  a  stratified 
wake  similar  to  that  expected  behind  a  submarine. 
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